/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */

//
// 文件: GPUConcreteKernels_impl.cuh
// 创建时间: 2025年10月4日
// 说明: GPU钢筋混凝土本构核函数实现（内联实现）
// 作者: OpenSees GPU加速项目组
//
// 注意：此文件由GPUConcreteKernels.h包含，提供设备函数的inline实现
// 算法参考:
//   - Lee & Fenves (1998): Plastic-Damage Model for Cyclic Loading of Concrete
//   - Lubliner et al. (1989): A Plastic-Damage Model for Concrete
//

// 不要在这里include GPUConcreteKernels.h，避免循环包含

//==============================================================================
// GPU辅助设备函数实现
//==============================================================================

/**
 * @brief 计算弹性刚度矩阵（3D各向同性）
 */
__device__ __inline__ void computeElasticStiffness3D(
    double E,
    double nu,
    double* Ce
)
{
    // 材料常数
    double G = E / (2.0 * (1.0 + nu));      // 剪切模量
    double K = E / (3.0 * (1.0 - 2.0 * nu)); // 体积模量

    // 初始化为零
    for (int i = 0; i < 36; i++) {
        Ce[i] = 0.0;
    }

    // Voigt符号的弹性刚度矩阵（6x6）
    // Ce = K * Ivp + 2G * Id
    // 其中 Ivp(i,j) = 1 (i,j < 3), Id 为偏应力投影

    // Ivp贡献（体积部分）
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            Ce[i * 6 + j] += K;
        }
    }

    // Id贡献（偏应力部分）
    for (int i = 0; i < 6; i++) {
        if (i < 3) {
            Ce[i * 6 + i] += 2.0 * G;
            for (int j = 0; j < 3; j++) {
                Ce[i * 6 + j] -= 2.0 * G / 3.0;
            }
        } else {
            Ce[i * 6 + i] += G; // 剪切分量
        }
    }
}

/**
 * @brief 应力分解函数（代数法）
 */
__device__ __inline__ void stressDecomposition(
    const double* stress,
    double* stress_pos,
    double* stress_neg,
    double* Q_pos,
    double* Q_neg
)
{
    const double tol = 1e-8;

    // 初始化投影矩阵为零
    for (int i = 0; i < 36; i++) {
        Q_pos[i] = 0.0;
        Q_neg[i] = 0.0;
    }

    // 代数分解（Voigt符号）
    for (int i = 0; i < 6; i++) {
        if (stress[i] > tol) {
            stress_pos[i] = stress[i];
            stress_neg[i] = 0.0;
            Q_pos[i * 6 + i] = 1.0;
            Q_neg[i * 6 + i] = 0.0;
        } else if (stress[i] < -tol) {
            stress_pos[i] = 0.0;
            stress_neg[i] = stress[i];
            Q_pos[i * 6 + i] = 0.0;
            Q_neg[i * 6 + i] = 1.0;
        } else {
            stress_pos[i] = stress[i] / 2.0;
            stress_neg[i] = stress[i] / 2.0;
            Q_pos[i * 6 + i] = 0.5;
            Q_neg[i * 6 + i] = 0.5;
        }
    }
}

/**
 * @brief 计算应力不变量（八面体应力）
 */
__device__ __inline__ void computeStressInvariants(
    const double* stress,
    double& sigoct,
    double& tauoct
)
{
    // 八面体正应力
    sigoct = (stress[0] + stress[1] + stress[2]) / 3.0;

    // 第二偏应力不变量J2
    double J2 = ((stress[0] - stress[1]) * (stress[0] - stress[1]) +
                 (stress[0] - stress[2]) * (stress[0] - stress[2]) +
                 (stress[1] - stress[2]) * (stress[1] - stress[2])) / 6.0 +
                stress[3] * stress[3] + stress[4] * stress[4] + stress[5] * stress[5];

    // 八面体剪应力
    tauoct = sqrt(2.0 * J2 / 3.0);
}

/**
 * @brief 计算应力范数
 */
__device__ double computeStressNorm(const double* stress)
{
    double norm = sqrt(stress[0] * stress[0] +
                       stress[1] * stress[1] +
                       stress[2] * stress[2] +
                       2.0 * stress[3] * stress[3] +
                       2.0 * stress[4] * stress[4] +
                       2.0 * stress[5] * stress[5]);
    return norm;
}

/**
 * @brief 矩阵向量乘法（6x6矩阵 × 6向量）
 */
__device__ __inline__ void matrixVectorMult6x6(
    const double* matrix,
    const double* vector,
    double* result
)
{
    for (int i = 0; i < 6; i++) {
        result[i] = 0.0;
        for (int j = 0; j < 6; j++) {
            result[i] += matrix[i * 6 + j] * vector[j];
        }
    }
}

/**
 * @brief 向量点积（6维）
 */
__device__ double vectorDot6(const double* v1, const double* v2)
{
    double dot = 0.0;
    for (int i = 0; i < 6; i++) {
        dot += v1[i] * v2[i];
    }
    return dot;
}

/**
 * @brief 张量积运算（向量×向量 → 矩阵）
 */
__device__ __inline__ void tensorProduct6x6(
    const double* v1,
    const double* v2,
    double* result
)
{
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            result[i * 6 + j] = v1[i] * v2[j];
        }
    }
}

/**
 * @brief 矩阵求逆（简化版：仅对角阵求逆）
 */
__device__ int matrixInverse6x6(double* matrix)
{
    // 简化实现：假设对角矩阵
    const double tol = 1e-14;

    for (int i = 0; i < 6; i++) {
        double diag = matrix[i * 6 + i];
        if (fabs(diag) < tol) {
            return -1; // 奇异矩阵
        }
        matrix[i * 6 + i] = 1.0 / diag;
    }

    return 0;
}

//==============================================================================
// GPU塑性损伤混凝土本构实现
//==============================================================================

/**
 * @brief GPU塑性损伤混凝土应力更新核函数
 *
 * 注意：使用__inline__以支持跨编译单元调用
 */
__device__ __inline__ void computePlasticDamageConcrete3D(
    const double* strain_increment,
    const double* strain_total,
    const PlasticDamageConcreteParams& params,
    ConcreteState& state,
    double* stress,
    double* D_tangent
)
{
    const double tol = 1.0e-5;

    // 提取材料参数
    double E = params.E;
    double nu = params.nu;
    double ft = params.ft;
    double fc = params.fc;
    double beta = params.beta;
    double Ap = params.Ap;
    double An = params.An;
    double Bn = params.Bn;

    // 计算参数k和初始损伤阈值
    double f2c = 1.16 * fc;
    double k = sqrt(2.0) * (f2c - fc) / (2.0 * f2c - fc);
    double rp0 = ft / sqrt(E);
    double rn0 = sqrt((-k + sqrt(2.0)) * fc / sqrt(3.0));

    // 计算弹性刚度矩阵Ce
    double Ce[36];
    computeElasticStiffness3D(E, nu, Ce);

    // ====== 第一部分：塑性计算 ======

    // 计算试探有效应力
    double Ce_Deps[6];
    matrixVectorMult6x6(Ce, strain_increment, Ce_Deps);

    double sige_tr[6];
    for (int i = 0; i < 6; i++) {
        sige_tr[i] = state.sige[i] + Ce_Deps[i];
    }

    // 应力分解
    double sigpos[6], signeg[6];
    double Qpos[36], Qneg[36];
    stressDecomposition(sige_tr, sigpos, signeg, Qpos, Qneg);

    // 计算负应力不变量
    double sigoct, tauoct;
    computeStressInvariants(signeg, sigoct, tauoct);
    double taun = sqrt(sqrt(3.0) * (k * sigoct + tauoct));

    // 塑性检查和修正
    double Cbar[36];
    double sige[6];

    if ((taun - state.rn) <= (tol * rn0)) {
        // 弹性状态
        for (int i = 0; i < 6; i++) {
            sige[i] = sige_tr[i];
        }
        for (int i = 0; i < 36; i++) {
            Cbar[i] = Ce[i];
        }
    } else {
        // 塑性状态 - Return Mapping
        double nrm = computeStressNorm(sige_tr);
        if (nrm < 1e-14) {
            for (int i = 0; i < 6; i++) sige[i] = 0.0;
            for (int i = 0; i < 36; i++) Cbar[i] = Ce[i];
        } else {
            double L_tr[6];
            for (int i = 0; i < 6; i++) {
                L_tr[i] = sige_tr[i] / nrm;
            }

            double L_trDotDeps = vectorDot6(L_tr, strain_increment);

            // 更新塑性应变
            for (int i = 0; i < 6; i++) {
                double Depse_tr = strain_total[i] - state.eps_p[i];
                state.eps_p[i] += beta * E * L_trDotDeps / nrm * Depse_tr;
            }

            double lam = 1.0 - beta * E / nrm * L_trDotDeps;

            // 修正的有效应力
            for (int i = 0; i < 6; i++) {
                sige[i] = lam * sige_tr[i];
            }

            // 重新检查损伤
            stressDecomposition(sige, sigpos, signeg, Qpos, Qneg);
            computeStressInvariants(signeg, sigoct, tauoct);
            taun = sqrt(sqrt(3.0) * (k * sigoct + tauoct));

            if ((taun - state.rn <= tol * rn0) || (L_trDotDeps <= 0)) {
                for (int i = 0; i < 6; i++) {
                    sige[i] = sige_tr[i];
                }
                for (int i = 0; i < 36; i++) {
                    Cbar[i] = Ce[i];
                }
            } else {
                // 计算切线刚度（简化）
                for (int i = 0; i < 36; i++) {
                    Cbar[i] = lam * Ce[i];
                }
            }
        }
    }

    // 更新有效应力状态
    for (int i = 0; i < 6; i++) {
        state.sige[i] = sige[i];
    }

    // ====== 第二部分：损伤计算 ======

    // 最终应力分解
    stressDecomposition(sige, sigpos, signeg, Qpos, Qneg);

    // 计算等效应力 taup（拉伸）
    double Ce_inv_sigpos[6];
    for (int i = 0; i < 6; i++) {
        Ce_inv_sigpos[i] = sigpos[i] / E; // 简化：假设各向同性
    }
    double taup = sqrt(vectorDot6(sigpos, Ce_inv_sigpos));

    // 计算等效应力 taun（压缩）
    computeStressInvariants(signeg, sigoct, tauoct);
    taun = sqrt(sqrt(3.0) * (k * sigoct + tauoct));

    // 拉伸损伤
    double Ddp_Drp = 0.0;
    if ((taup - state.rp) > (tol * rp0)) {
        state.rp = taup;
        state.dp = 1.0 - rp0 / state.rp * exp(Ap * (1.0 - state.rp / rp0));
        Ddp_Drp = (Ap * state.rp + rp0) / (state.rp * state.rp) *
                  exp(Ap * (1.0 - state.rp / rp0));

        state.dp *= (1.0 - tol);
        Ddp_Drp *= (1.0 - tol);

        if (state.dp > 1.0 - tol) {
            state.dp = 1.0 - tol;
            Ddp_Drp = 0.0;
        }
    }

    // 压缩损伤
    double Ddn_Drn = 0.0;
    if ((taun - state.rn) > (tol * rn0)) {
        state.rn = taun;
        state.dn = 1.0 - rn0 / state.rn * (1.0 - An) -
                   An * exp(Bn * (1.0 - state.rn / rn0));
        Ddn_Drn = rn0 / (state.rn * state.rn) * (1.0 - An) +
                  An * Bn / rn0 * exp(Bn * (1.0 - state.rn / rn0));

        state.dn *= (1.0 - tol);
        Ddn_Drn *= (1.0 - tol);

        if (state.dn > 1.0 - tol) {
            state.dn = 1.0 - tol;
            Ddn_Drn = 0.0;
        }
    }

    // ====== 第三部分：应力更新和切线刚度 ======

    // 最终应力
    for (int i = 0; i < 6; i++) {
        stress[i] = (1.0 - state.dp) * sigpos[i] + (1.0 - state.dn) * signeg[i];
    }

    // 切线刚度矩阵（简化计算）
    // C = (1-dp)*Qpos*Cbar + (1-dn)*Qneg*Cbar
    double Dsigpos_Deps[36], Dsigneg_Deps[36];

    // Dsigpos_Deps = Qpos * Cbar
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            double sum = 0.0;
            for (int k = 0; k < 6; k++) {
                sum += Qpos[i * 6 + k] * Cbar[k * 6 + j];
            }
            Dsigpos_Deps[i * 6 + j] = sum;
        }
    }

    // Dsigneg_Deps = Qneg * Cbar
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < 6; j++) {
            double sum = 0.0;
            for (int k = 0; k < 6; k++) {
                sum += Qneg[i * 6 + k] * Cbar[k * 6 + j];
            }
            Dsigneg_Deps[i * 6 + j] = sum;
        }
    }

    // 最终切线刚度（简化：忽略损伤梯度项）
    for (int i = 0; i < 36; i++) {
        D_tangent[i] = (1.0 - state.dp) * Dsigpos_Deps[i] +
                       (1.0 - state.dn) * Dsigneg_Deps[i];
    }
}

/**
 * @brief FAReinforcedConcrete应力更新（简化实现）
 */
__device__ __inline__ void computeFAReinforcedConcrete(
    const double* strain_increment,
    const FAReinforcedConcreteParams& params,
    ConcreteState& state,
    double* stress,
    double* D_tangent
)
{
    // TODO: 实现FAReinforcedConcrete材料
    // 当前使用弹性模型作为占位符
    double Ce[36];
    computeElasticStiffness3D(params.E, params.nu, Ce);

    matrixVectorMult6x6(Ce, strain_increment, stress);

    for (int i = 0; i < 36; i++) {
        D_tangent[i] = Ce[i];
    }
}

/**
 * @brief ConcreteL01应力更新（简化实现）
 */
__device__ __inline__ void computeConcreteL01(
    const double* strain_increment,
    const ConcreteL01Params& params,
    ConcreteState& state,
    double* stress,
    double* D_tangent
)
{
    // TODO: 实现ConcreteL01材料
    // 当前使用弹性模型作为占位符
    double nu = 0.2; // 假设泊松比
    double Ce[36];
    computeElasticStiffness3D(params.E, nu, Ce);

    matrixVectorMult6x6(Ce, strain_increment, stress);

    for (int i = 0; i < 36; i++) {
        D_tangent[i] = Ce[i];
    }
}
